{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Volatility Forecasting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_This setup code is required to run in an IPython notebook_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "sns.set_style(\"darkgrid\")\n",
    "plt.rc(\"figure\", figsize=(16, 6))\n",
    "plt.rc(\"savefig\", dpi=90)\n",
    "plt.rc(\"font\", family=\"sans-serif\")\n",
    "plt.rc(\"font\", size=14)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data\n",
    "These examples make use of S&P 500 data from Yahoo! that is available from `arch.data.sp500`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "import arch.data.sp500\n",
    "from arch import arch_model\n",
    "\n",
    "data = arch.data.sp500.load()\n",
    "market = data[\"Adj Close\"]\n",
    "returns = 100 * market.pct_change().dropna()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basic Forecasting\n",
    "Forecasts can be generated for standard GARCH(p,q) processes using any of the three forecast generation methods:\n",
    "\n",
    "* Analytical\n",
    "* Simulation-based\n",
    "* Bootstrap-based\n",
    "\n",
    "Be default forecasts will only be produced for the final observation in the sample so that they are out-of-sample.\n",
    "\n",
    "Forecasts start with specifying the model and estimating parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "am = arch_model(returns, vol=\"Garch\", p=1, o=0, q=1, dist=\"Normal\")\n",
    "res = am.fit(update_freq=5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "forecasts = res.forecast()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Forecasts are contained in an `ARCHModelForecast` object which has 4 attributes:\n",
    "\n",
    "* `mean` - The forecast means\n",
    "* `residual_variance` - The forecast residual variances, that is $E_t[\\epsilon_{t+h}^2]$\n",
    "* `variance` - The forecast variance of the process, $E_t[r_{t+h}^2]$.  The variance will differ from the residual variance whenever the model has mean dynamics, e.g., in an AR process.\n",
    "* `simulations` - An object that contains detailed information about the simulations used to generate forecasts.  Only used if the forecast `method` is set to `'simulation'` or `'bootstrap'`.  If using `'analytical'` (the default), this is `None`.\n",
    "\n",
    "The three main outputs are all returned in `DataFrame`s with columns of the form `h.#` where `#` is the number of steps ahead.  That is, `h.1` corresponds to one-step ahead forecasts while `h.10` corresponds to 10-steps ahead.\n",
    "\n",
    "The default forecast only produces 1-step ahead forecasts."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(forecasts.mean.iloc[-3:])\n",
    "print(forecasts.residual_variance.iloc[-3:])\n",
    "print(forecasts.variance.iloc[-3:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Longer horizon forecasts can be computed by passing the parameter `horizon`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "forecasts = res.forecast(horizon=5)\n",
    "print(forecasts.residual_variance.iloc[-3:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Alternative Forecast Generation Schemes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fixed Window Forecasting\n",
    "\n",
    "Fixed-windows forecasting uses data up to a specified date to generate all forecasts after that date. This can be implemented by passing the entire data in when initializing the model and then using ``last_obs`` when calling ``fit``.  ``forecast()`` will, by default, produce forecasts after this final date.\n",
    "\n",
    "**Note** ``last_obs`` follow Python sequence rules so that the actual date in ``last_obs`` is not in the sample."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res = am.fit(last_obs=\"2011-1-1\", update_freq=5)\n",
    "forecasts = res.forecast(horizon=5)\n",
    "print(forecasts.variance.dropna().head())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Rolling Window Forecasting\n",
    "\n",
    "Rolling window forecasts use a fixed sample length and then produce one-step from the final observation.  These can be implemented using ``first_obs`` and ``last_obs``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "index = returns.index\n",
    "start_loc = 0\n",
    "end_loc = np.where(index >= \"2010-1-1\")[0].min()\n",
    "forecasts = {}\n",
    "for i in range(20):\n",
    "    sys.stdout.write(\".\")\n",
    "    sys.stdout.flush()\n",
    "    res = am.fit(first_obs=i, last_obs=i + end_loc, disp=\"off\")\n",
    "    temp = res.forecast(horizon=3).variance\n",
    "    fcast = temp.iloc[0]\n",
    "    forecasts[fcast.name] = fcast\n",
    "print()\n",
    "print(pd.DataFrame(forecasts).T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Recursive Forecast Generation\n",
    "\n",
    "Recursive is similar to rolling except that the initial observation does not change.  This can be easily implemented by dropping the ``first_obs`` input."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "index = returns.index\n",
    "start_loc = 0\n",
    "end_loc = np.where(index >= \"2010-1-1\")[0].min()\n",
    "forecasts = {}\n",
    "for i in range(20):\n",
    "    sys.stdout.write(\".\")\n",
    "    sys.stdout.flush()\n",
    "    res = am.fit(last_obs=i + end_loc, disp=\"off\")\n",
    "    temp = res.forecast(horizon=3).variance\n",
    "    fcast = temp.iloc[0]\n",
    "    forecasts[fcast.name] = fcast\n",
    "print()\n",
    "print(pd.DataFrame(forecasts).T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## TARCH"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Analytical Forecasts\n",
    "\n",
    "All ARCH-type models have one-step analytical forecasts.  Longer horizons only have closed forms for specific models.  TARCH models do not have closed-form (analytical) forecasts for horizons larger than 1, and so simulation or bootstrapping is required.  Attempting to produce forecasts for horizons larger than 1 using `method='analytical'` results in a `ValueError`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TARCH specification\n",
    "am = arch_model(returns, vol=\"GARCH\", power=2.0, p=1, o=1, q=1)\n",
    "res = am.fit(update_freq=5)\n",
    "forecasts = res.forecast()\n",
    "print(forecasts.variance.iloc[-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simulation Forecasts\n",
    "\n",
    "When using simulation- or bootstrap-based forecasts, an additional attribute of an `ARCHModelForecast` object is meaningful -- `simulation`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "fig, ax = plt.subplots(1, 1)\n",
    "var_2016 = res.conditional_volatility[\"2016\"] ** 2.0\n",
    "subplot = var_2016.plot(ax=ax, title=\"Conditional Variance\")\n",
    "subplot.set_xlim(var_2016.index[0], var_2016.index[-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "forecasts = res.forecast(horizon=5, method=\"simulation\")\n",
    "sims = forecasts.simulations\n",
    "\n",
    "x = np.arange(1, 6)\n",
    "lines = plt.plot(x, sims.residual_variances[-1, ::5].T, color=\"#9cb2d6\", alpha=0.5)\n",
    "lines[0].set_label(\"Simulated path\")\n",
    "line = plt.plot(x, forecasts.variance.iloc[-1].values, color=\"#002868\")\n",
    "line[0].set_label(\"Expected variance\")\n",
    "plt.gca().set_xticks(x)\n",
    "plt.gca().set_xlim(1, 5)\n",
    "legend = plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bootstrap Forecasts\n",
    "\n",
    "Bootstrap-based forecasts are nearly identical to simulation-based forecasts except that the values used to simulate the process are computed from historical data rather than using the assumed distribution of the residuals. Forecasts produced using this method also return an `ARCHModelForecastSimulation` containing information about the simulated paths."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "forecasts = res.forecast(horizon=5, method=\"bootstrap\")\n",
    "sims = forecasts.simulations\n",
    "\n",
    "lines = plt.plot(x, sims.residual_variances[-1, ::5].T, color=\"#9cb2d6\", alpha=0.5)\n",
    "lines[0].set_label(\"Simulated path\")\n",
    "line = plt.plot(x, forecasts.variance.iloc[-1].values, color=\"#002868\")\n",
    "line[0].set_label(\"Expected variance\")\n",
    "plt.gca().set_xticks(x)\n",
    "plt.gca().set_xlim(1, 5)\n",
    "legend = plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Value-at-Risk Forecasting\n",
    "Value-at-Risk (VaR) forecasts from GARCH models depend on the conditional mean, the conditional volatility and the quantile of the standardized residuals,\n",
    "\n",
    "$$VaR_{t+1|t} = -\\mu_{t+1|t} - \\sigma_{t+1|t} q_{\\alpha}$$\n",
    "\n",
    "where $q_{\\alpha}$ is the $\\alpha$ quantile of the standardized residuals, e.g., 5%. \n",
    "\n",
    "The quantile can be either computed from the estimated model density or computed using the empirical distribution of the standardized residuals.  The example below shows both methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "am = arch_model(returns, vol=\"Garch\", p=1, o=0, q=1, dist=\"skewt\")\n",
    "res = am.fit(disp=\"off\", last_obs=\"2017-12-31\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parametric VaR\n",
    "First, we use the model to estimate the VaR.  The quantiles can be computed using the `ppf` method of the distribution attached to the model. The quantiles are printed below.\n",
    "\n",
    "**Note**: `forecast` is called with `align=\"target\"` so that the forecasts are already aligned with the target and so do not need further shifting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "forecasts = res.forecast(align=\"target\")\n",
    "cond_mean = forecasts.mean[\"2018\":].dropna()\n",
    "cond_var = forecasts.variance[\"2018\":].dropna()\n",
    "q = am.distribution.ppf([0.01, 0.05], res.params[-2:])\n",
    "print(q)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we plot the two VaRs along with the returns.  The returns that violate the VaR forecasts are highlighted. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "value_at_risk = -cond_mean.values - np.sqrt(cond_var).values * q[None, :]\n",
    "value_at_risk = pd.DataFrame(value_at_risk, columns=[\"1%\", \"5%\"], index=cond_var.index)\n",
    "ax = value_at_risk.plot(legend=False)\n",
    "xl = ax.set_xlim(value_at_risk.index[0], value_at_risk.index[-1])\n",
    "rets_2018 = returns[\"2018\":].copy()\n",
    "rets_2018.name = \"S&P 500 Return\"\n",
    "c = []\n",
    "for idx in value_at_risk.index:\n",
    "    if rets_2018[idx] > -value_at_risk.loc[idx, \"5%\"]:\n",
    "        c.append(\"#000000\")\n",
    "    elif rets_2018[idx] < -value_at_risk.loc[idx, \"1%\"]:\n",
    "        c.append(\"#BB0000\")\n",
    "    else:\n",
    "        c.append(\"#BB00BB\")\n",
    "c = np.array(c, dtype=\"object\")\n",
    "labels = {\n",
    "    \"#BB0000\": \"1% Exceedence\",\n",
    "    \"#BB00BB\": \"5% Exceedence\",\n",
    "    \"#000000\": \"No Exceedence\",\n",
    "}\n",
    "markers = {\"#BB0000\": \"x\", \"#BB00BB\": \"s\", \"#000000\": \"o\"}\n",
    "for color in np.unique(c):\n",
    "    sel = c == color\n",
    "    ax.scatter(\n",
    "        rets_2018.index[sel],\n",
    "        -rets_2018.loc[sel],\n",
    "        marker=markers[color],\n",
    "        c=c[sel],\n",
    "        label=labels[color],\n",
    "    )\n",
    "ax.set_title(\"Parametric VaR\")\n",
    "leg = ax.legend(frameon=False, ncol=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Filtered Historical Simulation\n",
    "\n",
    "Next, we use the empirical distribution of the standardized residuals to estimate the quantiles.  These values are very similar to those estimated using the assumed distribution.  The plot below is identical except for the slightly different quantiles."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "std_rets = (returns[:\"2017\"] - res.params[\"mu\"]) / res.conditional_volatility\n",
    "std_rets = std_rets.dropna()\n",
    "q = std_rets.quantile([0.01, 0.05])\n",
    "print(q)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "value_at_risk = -cond_mean.values - np.sqrt(cond_var).values * q.values[None, :]\n",
    "value_at_risk = pd.DataFrame(value_at_risk, columns=[\"1%\", \"5%\"], index=cond_var.index)\n",
    "ax = value_at_risk.plot(legend=False)\n",
    "xl = ax.set_xlim(value_at_risk.index[0], value_at_risk.index[-1])\n",
    "rets_2018 = returns[\"2018\":].copy()\n",
    "rets_2018.name = \"S&P 500 Return\"\n",
    "c = []\n",
    "for idx in value_at_risk.index:\n",
    "    if rets_2018[idx] > -value_at_risk.loc[idx, \"5%\"]:\n",
    "        c.append(\"#000000\")\n",
    "    elif rets_2018[idx] < -value_at_risk.loc[idx, \"1%\"]:\n",
    "        c.append(\"#BB0000\")\n",
    "    else:\n",
    "        c.append(\"#BB00BB\")\n",
    "c = np.array(c, dtype=\"object\")\n",
    "for color in np.unique(c):\n",
    "    sel = c == color\n",
    "    ax.scatter(\n",
    "        rets_2018.index[sel],\n",
    "        -rets_2018.loc[sel],\n",
    "        marker=markers[color],\n",
    "        c=c[sel],\n",
    "        label=labels[color],\n",
    "    )\n",
    "ax.set_title(\"Filtered Historical Simulation VaR\")\n",
    "leg = ax.legend(frameon=False, ncol=3)"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
